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Synchronization is an important and prevalent phenomenon in natural and engineered systems. 

In many dynamical networks, the coupling is balanced or adjusted in order to admit global synchro¬ 
nization, a condition called Laplacian coupling. Many networks exhibit incomplete synchronization, 
where two or more clusters of synchronization persist, and computational group theory has recently 
proved to be valuable in discovering these cluster states based upon the topology of the network. In 
the important case of Laplacian coupling, additional synchronization patterns can exist that would 
not be predicted from the group theory analysis alone. The understanding of how and when clus¬ 
ters form, merge, and persist is essential for understanding collective dynamics, synchronization, 
and failure mechanisms of complex networks such as electric power grids, distributed control net¬ 
works, and autonomous swarming vehicles. We describe here a method to find and analyze all of the 
possible cluster synchronization patterns in a Laplacian-coupled network, by applying methods of 
computational group theory to dynamically-equivalent networks. We present a general technique to 
evaluate the stability of each of the dynamically valid cluster synchronization patterns. Our results 
are validated in an electro-optic experiment on a 5 node network that confirms the synchronization 
patterns predicted by the theory. 
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I. INTRODUCTION 

Synchronization of oscillators in large networks has been an interesting problem for many years. It’s a phenomenon 
that shows up in many natural and man-made systems [1—3]. A particular type of synchronization is global syn¬ 
chronization where the oscillators all follow the same trajectory in state space. A network where this is desirable 
would be generators in a power grid, as well as some control networks and swarming autonomous vehicles. Although 
global synchronization has a well-developed theory [4-6] a more recently studied, more complex phenomenon we 
will call cluster synchronization, which we will abbreviate as CS, has attracted considerable attention [7-16]. In CS 
the network evolves into subsets of oscillators in which members of the same cluster are synchronized to the same 
trajectory, but members of different clusters are not. Such synchronized clusters may show up in swarms of animals, 
where the network is the simple visual link to one’s neighbors, or swarms of unmanned autonomous vehicles which 
are connected by a local communication network. Clusters may also show up in power grids where they would be a 
sign of a problem, that is, loss of global synchronization. 

Given the increase in man-made networks and the growing use of network theory to describe natural systems 
(e.g., food webs, neuronal and genetic networks) it is important to develop a basic approach to determining what 
cluster structures are possible in a given network. In this paper we show what methods can be used to do this using 
the concept of oscillators coupled through connection to other (not necessarily all) oscillators in the network. More 
importantly, we extend these methods to the currently unsolved problem of finding clusters in networks of oscillators 
which have a self-coupling to balance incoming signals from other oscillators, often called Laplacian coupling (more 
on this below). These allow synchronization clusters that elude other methods to find cluster patterns. We show how 
to use and extend symmetry methods to find all possible clusters in such networks of Laplacian-coupled oscillators. 
We first show how network symmetries can lead to cluster synchronization. We then show how we can go beyond this 
to analyze CS resulting from Laplacian Coupling. 

In Sec. II we review the concepts of symmetries and clusters in networks of coupled oscillators. In Sec. Ill we 
discuss methods to uncover all of the possible cluster synchronization patterns in a given network. Our main results 
are contained in Sec. IV, where we present a stability analysis that applies to any cluster synchronization pattern. 
We also present an elecro-optic experiment that confirms the patterns of synchronization predicted by the theory. 
Finally, the conclusions are presented in Sec. V. 



FIG. 1. (color online) (a) A network of 4 identical oscillators coupled through 3 identical links, (b) the same network after a 
reflection operation, (c) the same network after a rotation operation, (d) an 11-node network showing 3 clusters (blue, green, 
white). 
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II. SYMMETRIES AND CLUSTERS IN NETWORKS 

Fig. 1 (a) shows a 4-oscillator or 4-node network, where the oscillators are identical as are the couplings between 
them, which are bidirectional, meaning that signals flow in both directions to the oscillators by the same amount 
and influence on the connected oscillators. This network has a total of 6 symmetries. We show two. In Fig 1 (b) we 
show the result of a reflection (shown in (a)) interchanging nodes 1 and 2. The structure of the network remains the 
same. In Fig 1 (c) we show a rotation of the network by 120° which also leaves the network indistinguishable from 
the original in (a). These symmetries are manifest in the symmetries of a matrix that describes the network, the 
adjacency matrix. This matrix contains a 1 in each position of each row for each node that is connected to the node 
corresponding to the row number. For the network in Fig. 1 (a) the adjacency matrix A is, 


/0 0 0 1 \ 
[ 0 0 0 1 1 
0 0 0 1 
\1 1 1 0 / 


(1) 


The adjacency matrix plays a crucial role in modeling the dynamics of many networks of identical nodes, as it 
provides the coupling between the oscillators at each node. We will write the dynamics of the networks as follows, 


x; = F(x») + a'^2A ij H(x j ), (2) 

3 

i = 1, where x, is the vector of dynamical variables for the i th oscillator, x, is the time derivative of the ith 

node’s variables, N is the number of oscillators (nodes), F is the vector field of each node (governing each oscillator’s 
isolated dynamics), and H is a coupling function for each link in the network of of each oscillator to another. The 
entries Aij represent the strength of the coupling from node j to node i. Several papers [8, 11, 12, 15, 17-20] have 
used Eq. (2) to model the dynamics of a network. 

A network symmetry applied to the adjacency matrix leaves it unchanged. We recall that symmetries of an object 
(a network in this case) form a mathematical group Q. Any element in this group, say g is represented in the space of 
network nodes as a permutation matrix, say R g . The invariance of A under the action of the symmetry immediately 
implies R g A = AR g , A commutes with all group actions and the equations of motion for nodes that are mapped into 
each other are the same. For example, in the network in Fig. 1 (a), nodes 1,2, and 3 have the same equations of motion, 
so if they are started from the same initial conditions, they will remain synchronized indefinitely. Node 4 cannot be 
permuted into any of the other nodes and it will not synchronize with the others, hence we color it differently in the 
Figure parts (b) and (c). This is the intimate relationship between symmetry and dynamics in networks. We see that 
it immediately separates the network into two clusters 1,2,3 and 4. 

Many of the recent studies of cluster synchronization are on particular networks that are known or engineered to 
exhibit cluster synchrony. But as we showed above group theory provides a general approach to finding clusters in 
arbitrary networks. Steps toward more general approaches using group theory to analyze cluster synchronization 
began with the work of Golubitsky, Stewart, and Schaeffer [21, 22], where network symmetries are known and can be 
shown to support CS. Recently, computational methods have been used to study simple symmetric networks and an 
approach has been developed that relates the symmetries with the emergence of the CS states [23]. Finally, we showed 
that such approaches can be applied to more complex networks with hundreds of oscillators using computational group 
theory [24], 

An alternative description of the network dynamics is the following in the case of Laplacian coupling, 

X; = F(xj) ' H(x ; ) - H(xj)], (3) 

i 

where the coupling from oscillator j to oscillator i is given by the the difference between the output functions H(xj) 
and H(xj). Several papers [2, 5, 7, 10, 20, 23, 25] have used Eq. (3) to model the dynamics of a network. Equation 
(3) can be rewritten as follows, 

=F(x i ) + o-y^L ij H(xj), (4) 

3 

which has the same structure as Eq. (2) but now the adjacency matrix has been replaced by the Laplacian matrix, 
L = {Lij}, L ig = A^ — Sjj J2j where Sij is the Kronecker delta. By construction then we have that the sums of 
the rows of the matrix L are equal to zero, i.e. the inputs to the 7th node are balanced by the diagonal self-coupling. 
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In Ref. [26] Golubitsky, Stewart, and Torok have shown that for networks that have balanced coupling (all nodes 
receive the same cumulative input weights, accounting for adjacent nodes and self-coupling, an example of which is 
Eq. (4)), CS can emerge in many patterns that are not directly the result of symmetries. An example of this is global 
synchronization, which is not a result of symmetries in the network. 

In general, the patterns of cluster synchrony that can be observed in a network are not unique [27, 28], hence 
an important problem is that of determining the parameter ranges for stability and multi-stability for the observed 
patterns. While it is known that stability of the global synchronization state for an arbitrary network can be charac¬ 
terized by using the master stability function formalism [29] , a corresponding analysis that applies to the CS patterns 
that may emerge in a network is not available. In this paper we address this problem by providing necessary and 
sufficient conditions for stability of each individual CS pattern under very general assumptions. Our analysis applies 
to systems for which the functions describing the individual dynamics (possibly chaotic) and the interactions between 
the systems are arbitrary and to both the cases that the network is described by an adjacency matrix (Eq. (2)) or by 
a Laplacian matrix (Eq. (4)), both of which are used to model network interactions. 

In what follows we show how to find all of the CS patterns that may emerge in a given network topology, described 
either by an adjacency matrix or by a Laplacian matrix and how to evaluate the stability of each allowed pattern. 
We demonstrate all the above phenomena in an electro-optic experiment on a 5 node network that displays all of the 
possible CS patterns predicted by the theory. 


III. ANALYZING CLUSTER SYNCHRONIZATION PATTERNS 


Here we attempt at addressing the following problem: given a network structure (either in terms of an adjacency 
matrix in Eq. (2) or of a Laplacian matrix in Eq. (4)) can we find all of the cluster synchronization (CS) patterns 
that are allowed? For simplicity, we will proceed under the assumption that the network dynamics is described by 
Eq. (4) but all of our results include the simpler case that the dynamics is described by Eq. (2). Indeed, as we will 
see, the case of the Laplacian matrix is in general more complex to deal with than that of the adjacency matrix. So 
we consider the most difficult case. 


First of all, we note that a symmetry of the adjacency matrix is also a symmetry of the corresponding Laplacian 
matrix and viceversa [24]. Suppose Q is a group of permutations of the nodes of the network which leaves the coupling 
matrix L invariant. Then for each g £ Q we have a permutation matrix R g that operates on the set of all node vectors 
x = (x!..., xtv) T - Since R g L = LR g , this means d^RgX^/dt = F(R g x. i ) + <rJY LijH(R g x.j), i.e., the symmetry 
operation leaves the equations of motion unchanged. Hence, the subset of nodes permuted among each other by the 
group will remain synchronized if started in a synchronized state. We will refer to these subsets of nodes as clusters. 
The synchronized states for each cluster are flow invariant. 

An approach to the construction of all allowed CS patterns in a network has been proposed in [30]. While the 
method presented in [30] is general, as it applies to any network topology, it is computationally expensive. Here 
we argue that for the case of symmetric networks, a faster approach can often be followed that takes advantage of 
computational group theory, which is quite efficient. At each step we show the results of this method applied to a 
particular case. 


We first decompose the symmetry group Q into subgroups Hi, i = 1,..., v , each of which acts only on some subset 
of clusters (often only one), but not on any of the others [31, 32]. For this reason, we will refer to these subgroups as 
cluster groups and the original group is a direct product of all cluster groups Q = Hi x ... x H u . We further decompose 
each cluster group Hi into all of its possible subgroups, which will give us a natural set of symmetry-breaking paths. 
These subgroups provide the full range of possible symmetry clusters from the original full symmetry clusters to 
subclusters to the trivial case where each node is in its own cluster, i.e. no symmetries. While for the case of the 
adjacency matrix, this allows one to find all of the possible CS patterns, in the case of the Laplacian matrix, these 
patterns are certainly valid but there may be other valid patterns that are not predicted by the computational group 
theory analysis. Extra steps are thus required in order to find these additional patterns. 

We examine a particular case, a 5-node network (Fig. 2), with adjacency matrix 


A = 


/0 1 0 1 1 \ 
1 0 1 0 1 
0 10 11 
10 10 1 
Vi i i i 0 / 


( 5 ) 
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and Laplacian matrix 


^—3 10 11 


L = 
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FIG. 2. (color online) On the left: all possible patterns displayed when the network connectivity is given by the adjacency 
matrix (5). On the right: additional patterns displayed when the network connectivity is given by the Laplacian matrix (6). 

Figure 2 shows all the allowed patterns, where nodes belonging to the same cluster are colored the same. 

Cluster patterns Al to A5 in Fig. 2 can be found by performing an analysis of the group Q and all of its subgroups. 
These are all the patterns that may emerge for the adjacency matrix (5) and the only patterns that may emerge by 
symmetry for the Laplacian matrix (6). The orbits of the original symmetry group are {1,2,3,4} and {5} by itself 
and are associated with two cluster groups: Hi which permutes the first and H 2 which is only the identity for the 
2nd single node cluster. Other patterns (A2-A5) are possible based on subgroups of the original group. These are the 
result of symmetry breaking bifurcations. 

Now we create new potential clusters by first, choosing a set of cluster groups (Hi) and/or their subgroups (one 
for each cluster group). Together these determine a subgroup Q' of the original group Q. Second, we combine or 
merge some of these clusters as candidates for new synchronized clusters (it’s our choice which to try). Because we 
are merging clusters or subclusters from different original symmetry clusters, the resulting CS patterns will not be the 
result of symmetries but may be dynamically valid when the coupling is Laplacian. Third, we can set those dynamical 
variables Xj equal to others in the merged cluster and see if their equations of motion are the same (guaranteeing flow 
invariance). However, examining equations of motion by eye for large networks can become prohibitive. 

There is a more direct way of doing step three above using the power and efficiency of group theory and com¬ 
putational discrete algebra tools [33, 34]. Examining the coupling term in Eq. (4) when clusters are synchronized 
the diagonal feedback term (H(xj)) of a node will cancel the coupling terms of nodes from the same merged clus¬ 
ter. Hence, in the synchronized state the dynamics behaves as though nodes from the same merged cluster are not 
connected. We therefore form a dynamically equivalent coupling matrix which is the original Laplacian matrix with 
off-diagonal components between nodes from merged clusters set to 0 and the diagonal values then set to the negative 
of the new row sums. We then perform the cluster group decompositions and subgroup constructions on the new, 
dynamically equivalent Laplacian (note that if the original Laplacian matrix is symmetric, so are also the dynamically 
equivalent matrices obtained through this construction). If some set of subgroups of this new coupling matrix has 
symmetries yielding clusters which are our merged clusters, then their dynamics is flow invariant in the synchronized 
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state and the cluster merging is possible. In this case we call the new clusters Laplacian clusters. All the above 
can be automated in software. Hence, in Laplacian networks, even when a CS pattern is not the direct result of a 
symmetry of the original coupling matrix, it is the result of a symmetry of a dynamically equivalent coupling matrix. 
This is particularly relevant in terms of computational complexity. In fact, while the problem of finding all of the 
symmetry operations of a given matrix has not been proven to be polynomial, efficient discrete algebra routines have 
been developed that make these computations possible even for very large networks (e.g., the Internet at the AS level, 
for which N = 22332, see Ref. [32] ). This is in contrast with the method proposed in [30], not based on evaluation 
of the symmetries, which has been shown to become inefficient even for networks of moderate size, e.g., A = 15. 

In Fig. 2 the only possibility of merging is for node 5 to be absorbed into one of the remaining clusters in patterns 
Ai,A 2 , ..., A$, which is shown on the right hand side. Analyzing our dynamically equivalent matrix as above, we find 
three additional allowed CS patterns, L i, L 3 , and L 4 . These are not predicted by the original symmetry analysis but 
can emerge when the coupling matrix is Laplacian. Note that the potential L 2 patterns from merging of the center 
with one of the other clusters in A2 will not lead to a new synchronized cluster (the center has two inputs from the 
unsynchronized cluster while the corners have only one). The group analysis of the dynamically equivalent network 
agrees with this. 

To conclude, when the network is described by an adjacency matrix, a full characterization of all the dynamically 
valid patterns can be obtained by taking advantage of available computational discrete algebra tools [33, 34]. These 
output the cluster groups, "Hi, I-L 2 and a decomposition of each cluster group 77; into all of its possible 
subgroups, which provides a natural set of symmetry-breaking paths. As discussed above, available computational 
group theory routines perform these tasks very efficiently, when compared to other possible methods (e.g., [30]). For 
this case our approach is always better than the state of the art. When the coupling is in Laplacian form, all of the 
(symmetry-related) CS patterns of the associated adjacency matrix are mantained. Moreover, it is possible for some 
of the cluster groups 77i, 7 ^ 2 ,7A an d some of their subgroups to merge to form new dynamically valid clusters. In 
order to find all of the dynamically valid mergings, it can be helpful to use the clusters and subclusters provided by the 
symmetry analysis as the building blocks of our algorithm. If the number of these different clusters and subclusters is 
equal to fi, the number of tests that need to be performed (for pairs, triplets, and so on) is upper bounded by 7? M , the 
/ 4 -th Bell number, hence it grows combinatorially with /r. For all the networks that do not display many symmetries, 
i.e., for which /t A, our approach based on computational group theory will be much faster than the one presented 
in [30]. Hence, in general it is a good idea to preliminarily run an analysis of the symmetries of the network to assess 
whether /4 < A and choose the most convenient approach based on this outcome. It should be noted that there is no 
need for testing for mergings between subgroups of the same group, which reduces the complexity of the calculations. 
In general, finding all of the dynamically valid patterns for Laplacian networks may be substantially harder than for 
networks described by a symmetric adjacency matrix. Similar limitations were observed in [30]. 


IV. STABILITY ANALYSIS AND EXPERIMENTAL VALIDATION 

Stability of cluster synchronization has been investigated for phase oscillators [23] and Stuart-Landau oscillators 
[35] , and for lattices of coupled systems [27, 28] . However, a general approach to analyze stability and multistability of 
CS patterns in arbitrary networks has not been developed. In [24] we studied a particular CS pattern corresponding 
to a minimum number of clusters (i.e., maximal symmetry) for the case that the connectivity of the network is in the 
form of an adjacency matrix. 

We now develop variational equations for the merged-cluster system so we can calculate the stability of each one 
of the allowed CS patterns. While a number of papers [23, 26-28, 35-37] have dealt with cluster synchronization 
in networks, only Refs. [23, 35] have considered the problem of stability for particular dynamics of the individual 
systems. Refs. [27, 28, 37] have emphasized that for arbitrary systems this is a difficult problem. In order to analyze 
stability, we start with the subgroup Q' of the original group which generated the clusters we want to merge. It is 
formed by a direct product of the subgroups we choose to use in our merged system. Using C m to represent each 
cluster of nodes, m = 1,..., TV, where A=number of clusters in Q ', we have the variational equation of Eq. (2), 


<5x(f) 


" M 

Y E^ m) ® DF(s m (t)) 

_m= 1 
M 

+<J Y( LEim) ) ® DU ( S ™{t)) 

m —1 


<5x(f), 


( 7 ) 


where the An-dimensional vector 5x(f) 


[5xi(t) T , Sx 2 (t ) T ,..., (5xjv(t) T ] T , L= the Laplacian-coupling matrix, and 
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Pattern 

Quotient Network Dynamics (right hand side is mod27r) 

Transverse Perturbations 

Al: 

X\ = X2 = X3 = #4 

x [ +1 = [08 - a)X(x\) + oX(x\) + 5] 
xl +1 = [(/?- 4a)X(x t 5 ) + 4aX(x{) + <5] 

i f +1 = \P + a A] DX(x\)rf, A = -3, -5 

A2: 

Xl = X 2 & X3 = X4 

x i rl = [{P ~ 2 a)X(x\) + aX(x\) + aX(x%) + 5] 
x 2 +1 = [($ — 3a)X(xl) + 3 aX(x\) + <5] 

4 +1 = rC/9 — 3a)X(xl) + 3aX(x\) + 5} 

V i +1 = [0 - 4a] DX(x\)rA - *DI(x\)r& 
V 3 +1 = -crDX^xX)^! + [P- 4a] a-DX^I)^ 

A3: 

Xl = X3 & X 2 — X4 

*i +1 = \{P ~ 2a)X{x\) + <jX(x 2 ) + aX(x\) + 5] 
x 2 +1 = [(P ~ 2 a)X{x t 2 ) + aX(x *) + aX(x%) + <5] 
xl +1 = f(/3 — 4a)X(xl) + 2aX(x\) + 2 aX(x t 2 ) + 5] 

rf +1 = \J3 + A a] DX(x \)rf , A = -3 

A4: xi = X 3 

x = [(/3 — 3cr)X(xi) + gX{x 2 ) + oX{x 4) + gX{x 5 ) + 
xl +1 = [(p - 3 g)X[x\) + 2gX(x\) + gX(x\) + S] 
x 4 +1 — [(P ~ 3cr)X(x|) + 2aX(x\) + aX{x\) + S] 

= [(/? — 4 cr)X(fE 5 ) + 2aX(x\) + < 7 X^X 2 ) + crX{x\ ) + $] 

r ] t+1 = \P + Aa] DX(x\)rf, A = -3 

LI: xi = X 2 = ... = X 5 

x [ +1 = \px(x\) + 6 } 

rf +1 = [P + a A] DX(x\)r]\ A = -3, -5 

L3: 

Xl = X3 = X 5 &ZX 2 = X 4 

x\ +1 = [(/? - 2a)X(x\) + 2aX(x t 2 ) + < 5 ] 
x 2 +1 = \{P ~ 3a)X(x t 2 ) + 3 aX(x\) + 5 ] 

^ +1 = [p + a A] DX{ xi)r)L A = -3, -5 

77* +1 = [P + aA] DI{xi)r&, A = -3 

L4: xi = X 3 — X 5 

*i +1 = [{P ~ 2<j)X(x\) + uX(x\) + crX(x\) + 5] 
x^ +1 = [{P — 3a)X(x2) + 3aX(x\) + 5] 

*4 +1 = \{P — 3a)X(x t 4) + 3 aX(x\) + 5] 

T) t+1 = \p + a A] DX(x Drf, A = -3, -5 


TABLE I. Equations used to evaluate stability of all the allowed CS patterns for the system described by Eq. (10) and coupling 
matrix (6). 


£( m ) is an TV-dimensional diagonal indicator matrix for each cluster such that E is equal to 1 if i £ C m and is equal 
to 0, otherwise, i = 1,TV. Note that we must use the original Laplacian matrix and not the dynamically equivalent 
one, which is used only for detecting synchronization flow invariance. 

As we showed in [24] we can first block diagonalize the coupling matrix L using the irreducible representations 
(IRR) of Q', which yields the transformed coupling matrix L' for A3 of Fig. 2, 


/ 

-4 

-V2 


0 

0 ^ 



-V2 

-3 

-2 

0 

0 



V2 

-2 

-3 

0 

0 

(8) 


0 

0 

0 

—3 

0 


V 

0 

0 

0 

0 




The upper-left block represents the variations on the synchronization manifold. It is 3 x 3 because there are three 
different trajectories in A3 (three clusters). The lower-right block represents the variations in the transverse manifold. 
These are the perturbations that take the system out of synchrony so it is these whose stability we want to calculate. 
Suppose we merge the center node (5) with nodes 1 and 3 to form the first merged state shown in L3. Geometrically, 
what must happen is that the dimension of the synchronization manifold must decrease by one (from 3 to 2) and the 
transverse manifold must increase by 1 (from 2 to 3). 

To obtain the new coordinates on the synchronization manifold, we note that basis vectors of a cluster on the 
synchronization manifold have a 1 in the position of each node that belongs to the cluster. For example, in A3 
the cluster (1,3) will have a (unit) vector of the form (1,0,1,0,0) and (5) will have the (unit) vector of the form 
(0,0,0,0,1). The merged cluster (1,3,5) will have the synchronization manifold vector which is their sum, (1,0,1,0,1). 
The transformation of this new synchronization vector to the IRR coordinates of L' provides the new synchronization 
direction and its orthogonal complement provides the new transverse direction. We use these two new vectors to 
transform the 2x2 sub-block associated with the (1,3) and (5) clusters in the 3x3 synchronization block to reduce 
the 3x3 synchronization manifold and increase the transverse manifold. This results in the final variational equation 
for the L3 case, 


M 

[ j(m) ® DF ( S ™) + <tL"jW 0 Z?H(s m ) 

















new transverse direction 


ncbr° nizat l° n 


FIG. 3. Reduction of the dimension of the three-dimensional synchronization manifold. This shows schematically how the 
merging of clusters (1,3) and (5) produces a new synchronization direction in the (1,3) and (5) plane of the synchronization 
manifold along with a new transverse direction orthogonal to the new synchronization direction. 


where we have linearized about the new synchronized merged cluster states {s l7 f] is the vector of variations 

of all nodes transformed to the merged coordinates as above, DF and DH are the Jacobians of the nodes’ vector field 
and coupling function, respectively, are the transformed E( m ' ) and, 


L" 


/—3 V6 0 0 0 \ 

V6 -2 0 0 0 

0 0 -5 0 0 

0 0 0 -3 0 

V 0 0 0 0 -3/ 


(9) 


In Eq.(9) the new synchronization block (in the upper-left-hand corner) represents the new clusters (1,3,5) and (2,4) 
and the new transverse direction is associated with the new diagonal value —5. Obviously, this can be generalized 
to more complex cluster mergings. An example of cluster merging is shown geometrically in Fig. 3. We remark that 
this method can be used to analyze stability of any dynamically valid CS pattern, for which knowledge of the block- 
diagonalized matrix V is available and it applies to both cases that the connectivity is given by an adjacency matrix 
or by a Laplacian matrix. The generalization of the above procedure simply uses the synchronization vectors for all 
the new clusters as the rows of a matrix for which all other components are zero. An application of a singular value 
decomposition then gives a basis for the synchronization block (the original synchronization vectors) and a basis for 
the orthogonal compliment which represents all the transverse directions. In some cases the latter can be simplified 
by a block diagonalization when all members of the block are in the same merged cluster. This makes it possible to 
automate this block-diagonalization to evaluate stability of all the CS patterns that can emerge in a given network 
topology. 

We show symmetry-breaking and the existence of Laplacian clusters using the experimental system described in 
detail in [24, 38] . This system employs a spatial light modulator (SLM) and a camera in a feedback configuration. The 
camera has a focal plane array (FPA) of 320 x 256 pixels, and an area of 8 x 6.4 mm 2 . The SLM has a resolution of 
512 x 512 pixels, and an active area of 7.68 x 7.68 mm 2 . A light emitting diode with a wavelength of 1550 nm is used 
to illuminate the modulator. The light passes through a polarizing beam-splitter, and a quarter-wave plate (QWP), 
so that circularly polarized light is incident on the SLM. The SLM imparts a programmable spatially-varying phase 
shift x between the two polarization components of the reflected light. The reflected light passes through the QWP 
and polarizer, and is imaged onto a 256 x 256 pixel square region of the camera’s FPA. The relationship between the 
phase shift x applied by the SLM and the normalized intensity I recorded by the camera is I(x) = (1 — cos a;) /2. 

Each oscillator corresponds to a square patch of 16 x 16 pixels on the SLM, which is imaged onto an 8 x 8 pixel 
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region of the camera’s FPA. Using a computer, the phase shift of the ith region, Xi, is iteratively updated based on 
the intensity measured by the camera according to: 


„t+1 


131 (xl) + a ^2 Lijl(xj) + S 


mod 2tt 


( 10 ) 


where cr is the coupling strength, which we will vary. Eq. (10) is a map version of Eq. (2) and can show fixed point, 
periodic or chaotic dynamics depending on the values of the parameters. Here, S = 0.525, /3 = 1.457T, and a was 
decreased from n to 0. 



FIG. 4. The figure shows the experimental synchronization error for each synchronization pattern as a function of the parameter 
cr. Underneath we plot the results of our stability analysis applied to each one of the cluster synchronization patterns, where 
a colored dot labels the values of sigma for which the corresponding pattern is stable. 


Fig. 4 shows the experimental synchronization error for each synchronization pattern shown in Fig. 2 as a function 
of a . The phases Xi were not reinitialized when a was updated. The synchronization error for each pattern is computed 
as (| Xi(t) — xf(t )|), where the symbol (•) indicates an average both in time and over the nodes * in a cluster and the 
clusters in a pattern, xf is the average state for the nodes in the cluster to which node i belongs. In the lower panel 
of Fig. 4, we plot the results of our stability analysis for each one of the CS patterns, where a colored dot labels the 
values of sigma for which the corresponding pattern is stable. In particular, a CS pattern was indicated to be stable 
when (i) all the numerically computed maximum Lyapunov exponents corresponding to the transverse blocks were 
found to be negative and (ii) the synchronous pattern was asymptotically valid, that is the CS pattern was observed 
after integrating its equations for a long time. The equations that were used to run these stability calculations are 
shown in Table I. 

In general, when two or more clusters merge into one, there are two independent effects on stability, as can be seen 
from the structure of the block-diagonalized matrix L": the first one is that the dimension of the synchronization 
block decreases, which determines the motion in the synchronization manifold; the second one is that new transverse 
blocks appear, with the other (pre-existing) transverse blocks remaining the same. As a consequence, since each one 
of these transverse blocks needs to be individually tested for stability, we expect that as the number of transversal 
blocks increases, the range of stability decreases. This is confirmed by our experimental results plotted in Fig. 4, 
showing that the a -range of stability becomes smaller for CS patterns that are characterized by higher symmetry. 
Exceptions to this rule are possible as the motion in the synchronization manifold (on which the transverse Lyapunov 
exponents depend) may also affect stability in ways that cannot be predicted by the analysis of the transverse blocks 
only. 
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FIG. 5. Experimental phase space plots with lines connecting successive iterates for a) three clusters (A3), b) two clusters (L3) 
and c) one cluster (LI). 


Fig. 5 shows the dynamics on the synchronization manifold for the symmetry pattern A3 and the two merged 
patterns, L3 and LI. In each transition A3 —> L3 —> LI the dimension of the synchronization manifold decreases by 
1 (from 3 to 2 to 1) and the transverse manifold increases by 1 (from 2 to 3 to 4). Fig. 6 shows three snapshots of 
the experimental dynamics for each one of the patterns A3,L3, and LI. Videos of the synchronization manifolds and 
patterns of Fig. 6 are available in the supplement. 



(a) 


(b) 


(c) 


FIG. 6. Snapshots of the experimental dynamics for the cases of a) three clusters (A3), b) two clusters (L3) and c) one cluster 
(LI). 


V. CONCLUSION 

We have studied the emergence of cluster synchronization in networks of coupled oscillators. First, we show how to 
obtain all of the possible dynamically valid CS patterns for an arbitrary network topology. The methods we illustrated 
above are general and will apply to nodes with any type of dynamics (ODEs, maps, etc.). This method also extends to 
directed networks, weighted networks, and labeled nodes (to represent different dynamics on each node), all of which 
are handled by the software package [34] . An important point is that these techniques work for any subgroups of the 
cluster groups which includes the trivial subgroups [24] (only the identity in each). Hence, we can approach merging 
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as a bottom up process and analyze any arbitrary merging of nodes into a cluster to determine if the dynamics allows 
a synchronized state and if it is stable in some parameter range. This means all possible clusters can be analyzed 
using our approach. We note that combining clusters as a “top down” approach would provide clusters that most 
likely are not from symmetries whereas bottom up would be a process that would include clusters easily obtained 
from symmetries, although it would be useful for cases where one has particular clusters in mind to analyze. 

Our main result is a technique to evaluate stability of all the dynamically valid CS patterns for both networks for 
which the connectivity is given by an adjacency matrix and by a Laplacian matrix. We predict that the range of 
stability typically becomes smaller for CS patterns that are characterized by higher symmetry, which is confirmed in 
our experimental system. 

We acknowledge help with computational group theory from Prof. David Joyner of the US Naval Academy and 
information and computer code from Ben D. MacArtlrur and Ruben J. Sanchez-Garcfa both of the University of 
Southampton. This work was funded by the Office of Naval Research. F.S. was supported by NSF grant CMMI- 
1400193. 
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